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Abstract 

We discuss the product of independent induced quaternion (/? = 4) Ginibre matrices, 
and the eigenvalue correlations of this product matrix. The joint probability density 
function for the eigenvalues of the product matrix is shown to be identical to that of a 
single Ginibre matrix, but with a more complicated weight function. We find the skew- 
orthogonal polynomials corresponding to the weight function of the product matrix, and 
use the method of skew-orthogonal polynomials to compute the eigenvalue correlation 
functions for product matrices of finite size. The radial behavior of the density of states 
is studied in the limit of large matrices, and the macroscopic density is discussed. The 
microscopic limit at the origin, at the edge(s) and in the bulk is also discussed for the 
radial behavior of the density of states. 
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1 Introduction 



Random matrix theory has found a vast number of applications in physics and mathematics as 
well as many other sciences; we refer to [1] for a recent compilation. But in certain applications 
it is insufficient to consider a single random matrix, and one must consider products of 
random matrices instead. Two limits are particular interesting: Products of a large number 
of matrices, and products of large matrices. As applications, which belong to the first category, 
we may mention: Lyapunov exponents for disordered and chaotic systems |2], multiplicative 
matrix- valued non-commutative diffusion processes [3], stability of ecological systems [3], and 
scattering of electromagnetic waves on random obstacles in telecommunication [5]. These 
have all been studied using products of random matrices. The limit of large matrices has 
been used to study: Transfer matrices in mesoscopic wires [6], large N c Wilson loops in 
Yang-Mills theory JTHSJE], and quantum chromo dynamics at finite chemical potential |1U] . 

Recently, several papers have appeared which discuss the product of large non-Hermitian 
matrices. The macroscopic, or mean, density of states for complex non-Hermitian matrices has 
been discussed using diagrammatic methods in 112], Proofs for the macroscopic behavior 
for products of complex non-Hermitian matrices can be found in |13[ 114] . The method of 
orthogonal polynomials was used in |15] to obtain explicit expressions for the eigenvalue 
correlation functions for a product of independent complex (/3 = 2) Ginibre matrices of finite 
size, and it was shown that the macroscopic limit agreed with previous results. A discussion 
of the microscopic behavior at the origin and at the edge as well as in the bulk of the spectrum 
was also presented in [15] , while hole probability and overcrowding at the origin was discussed 
in pi]. 

The product of two independent non-Hermitian matrices has been discussed for all three 
Ginibre ensembles: The real (/3 = 1) Ginibre ensemble was discussed in |17| IT8| IT9| 120] . the 
complex (/3 = 2) Ginibre ensemble in [21, 22j, and the quaternion (/3 = 4) Ginibre ensemble 
in [23]. But to our knowledge, the product of three or more Ginibre matrices has only been 
discussed in the complex case. 

In this paper we discuss the product of quaternion (/3 = 4) Ginibre matrices. The product 
of two quaternion Ginibre matrices is directly related to quantum chromodynamics with a 
finite chemical potential and fermions in the adjoint representation [23] . Here, we investigate 
the product of an arbitrary (finite) number of matrices. We will use the method of skew- 
orthogonal polynomials to obtain explicit expressions for the eigenvalue correlation functions 
for matrices of finite size. We will also discuss the radial behavior for the density of states 
in the limit of large matrices. This allows us to find the macroscopic density as well as 
microscopic limits at the origin, at the edge and in the bulk. 

We consider the product matrix 

= X1X2 ■ ■ ■ X n , (1) 

where Aj are independent matrices drawn from the induced quaternion (/? = 4) Ginibre 
ensemble, specified by the joint probability density function 

p£ df (A) = ^det(A) m exp[TrXtA], (2) 

where m is a non-negative number and Cj^ is a normalization constant. Note that we look at 
the induced Ginibre ensemble (if m 7^ 0); a similar generalization could be introduced in the 
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complex case [15J without much difficulty. Looking at the induced Ginibre ensemble allow us 
to introduce an inner edge for the spectrum, such that the eigenvalues are located within an 
annulus rather than a disk, as we will discuss more thoroughly in sections |4] and [5] The inner 
edge introduced by the parameter m was studied for a real (/? = 1) and a complex (/3 = 2) 
induced Ginibre matrix in |24j . That this situation is particular interesting is illustrated by 
the single ring theorem |25| [26] . 

Here and in the following X = [qij] denotes a non-Hermitian N x N matrix with indepen- 
dent quaternion entries 

q = g (0) + iq W + jg (2) + kg (3) , g (fc) G M. (3) 
Equivalently, we may represent the quaternion units as 2 x 2 matrices 

;«),,=(_»,;) «.*=(;;), 

where i is the imaginary unit. In this representation X = [qij] is a 2N x 2N matrix with a 
block structure given by 

/ u —v* 



V u 



u,veC. (5) 



Note that the determinant and the trace in the probability density function (J2J should be 
understood as applied to this 2N x 2N matrix. Moreover, it should be noted that the 
eigenvalues come in complex conjugate pairs such that the determinant in ([2]), and therefore 
the entire density function, is real and positive. For a more thorough discussion of quaternion 
matrices see e.g. [27J. 

The partition function corresponding to the product matrix ([I]) is given by 

/n . n 

|Z>X|I]p£df(*i) = / l^in^det^rexphTrX/Xi]. (6) 
i=i J i=i 

In our notation DX denotes the Euclidean volume element (e.g. the exterior product of all 
independent one-forms DX = [\ dX a ) and \DX\ is the corresponding unoriented volume 
element. This means that 

DX = f\DX i = f\ /\(dX,) ab = A A A (dX^) ab , (7) 

i=l i=l a,b i=l a,b q=0 

where % is the index for the n independent matrices, ab are the indices for the N 2 independent 
matrix entries, and q is the quaternion index. The corresponding unorientated volume element 
is 

n 3 

i^i= nnnc^w ^ 

i=l a,b q=0 

The main interest of this paper is twofold: First, we want to rewrite the partition function ([6]) 
in terms of the eigenvalues of product matrix P™, i.e. replace the joint probability density 
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functions, for the individual matrices, Xi, with joint probability density function, "P"'^ , 
for the eigenvalues of the product matrix, P™. Second, we want to use this insight to study 
the eigenvalue correlations of the product matrix. 

The paper is organized as follows: In section[2]we discuss how to obtain the joint probabil- 
ity density function, V^ d f, for the eigenvalues of the product matrix, P™. Section [3] explains 
how to find the eigenvalue correlations of the product matrix, P™, using the method of skew- 
orthogonal polynomials. The large- N behavior of the radial behavior for the density of states 
of the product matrix is discussed in sections [4] and [5] Section U] focuses on the macroscopic 
behavior, while section [5] discusses the microscopic limits. Some technical details are given in 
the appendices. 



2 Joint Probability Density Function 

In this section we outline the derivation of the joint probability density function, V™'^ , for the 
eigenvalues of the product matrix P™. We triangularize the independent Ginibre matrices, 
Xi, and use this parameterization to obtain the joint probability density function for the 
eigenvalues of the product matrix. Some technical details are given in appendix lAl 

We parameterize the quaternion matrices Aj (and therefore also product matrix P™) using 
a generalized Schur decomposition 

Xi = U i (A i + T l )Ur + \. (9) 

Here Aj are diagonal matrices, T, are strictly upper-triangular matrices and U{ are symplectic 
matrices with C/ n +i = U\. This parameterization was used to study the n = 2 case in [23], and 
a similar parameterization has been used to discuss the product of n independent complex 
(/? = 2) Ginibre matrices |15j . The product of two complex Ginibre matrices was first discussed 
in |21| 122] , Also the product of two independent real (/3 = 1) Ginibre matrices have been 
discussed in the literature |T71 [18j [20] . 

In the quaternion case, each Aj is a diagonal (block) matrix, where the diagonal entries 
are given by 

(Ai)aa = AL 0) + iAi? = jf), where x ia = xfj + iX^J E C; (10) 

each Tj is a strictly upper-triangular (block) matrix, e.g. (Tj) a fc = for a > b and (Ti) a b for 
a < b are given by the quaternion structure (0); and each J7j is a symplectic matrix. In order 
to ensure that the parameterization is unique it is necessary to introduce further constraints 
on the symplectic matrices, Ui. We see that the structure of each triangular matrix, Aj + Tj, 
is invariant under transformations 

Ui i y UiVi (11) 

where Vj is a diagonal matrix, where the diagonal elements are 2x2 diagonal unitary matrices. 
We choose Ui G Sp(N)/U(l) N such that the parameterization is unique. The degrees of 
freedom (d. o. f.) are (obviously) conserved by the parameterization, 

d.o.f.(Xj) = d.o.f.(A u T u U t ). (12) 

i=l,...,n 1=1, ...,n 
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It is well-known that for n = 1 the Jacobian corresponding to this change of variables is given 
by [28] 

N N 
a>b c=l 

where Zk are the eigenvalues of X% = P^L\- The derivation of the Jacobian for arbitrary n 
is given in appendix [X] The generalized Jacobian has a structure similar to the n = 1 case, 
except that the Jacobian now depends on the product of eigenvalues of the individual matrices 



Xi, i.e. J(n 



Moreover, it is straight forward to see that the joint probability density 



functions (j2|) for the individual matrices Xi are independent of the symplectic matrices, £7j, 
such that the integration over these degrees of freedom only contribute with a constant. Also 
the integration over the triangular matrices, Tj, contribute only with a constant; here the 
exponential term, exp[— TrT^Tj], ensures that the integral is finite. For this reason we may 
write the partition function ([6j in terms of the eigenvalues, Zk, of the product matrix ([T|, 



N 



k=l 



n / 



.zn) 



(14) 



and the joint probability density function for the product matrix is given by 



N 



Pjpdf (*!>••■ 



ZN) 



U<( Z c)\Zc 



N 
a>b 



\Zb 



(15) 



where C^J m is a normalization constant. We take the normalization to be C^ m 



1/4, which 

ensures that the density of states is normalized to the number of eigenvalues as we will see 
below. The weight function, w™(z), is defined using a Dirac delta function in order to write 
the joint probability function (|15[) as a function of the eigenvalues of the product matrix 
rather than a function of the eigenvalues of the individual matrices, Xi, 



n „ 

FT / d 2 Xi\x i \ 2m e-\^ 2 5 2 (z-x 1 X2---x n ) 



(16) 



For m = this weight function is identical to the weight discussed in [15], and here it was 
shown that the weight function could be expressed as a Meijer G-function. It turns out that 
the induced weight function (jl6|) also can be expressed as a Meijer G-function, 



<{z) 



7T 



0,n 



TO, 



, TO 



(17) 



The weight (j 1 7[) can be found following the same idea as in [15] . The main idea is to find a 
recursion relation for the Mellin transform of the family of weights {u>™} ra >i given by (j 1 6 1) . 
Solving the recursion relation and taking the inverse Mellin transform yields the weight func- 
tion given in equation (|17p . But since we already know the weight function for m = 0, we 
may take a shortcut. We use the fact that 



JV 



y[det(X t r = det(P™) = H\z k \ 2 , 



(18) 



i=l 



k=l 



5 



where Zk (and zt) are the eigenvalues of the product matrix P™ defined in <Q. It follows that 
the weight with m 7^ is related to weight with m = by 



<(*) = \z\ 2m w™=°(z). (19) 

It follows directly from the definition of the Meijer G-function |29j that the induced weight 
may be written as ([Tj 



3 Correlations and Skew- Orthogonal Polynomials 

One is often interested in the fc-point correlation functions for the eigenvalues, 

AT! N f 

Rk' m (zi,...,z k ) = JJ J^d 2 z h V^™(z 1 ,...,z N ), (20) 

which loosely speaking give the probability of finding k eigenvalues at the positions z\, . . . ,Zf. 
in the complex plane. In this section we will find an explicit expression for the eigenvalue 
correlations (J2UJ) using the method of skew-orthogonal polynomials. 

Following [30J, we use a standard identity for determinants to write the joint probability 
density function (J15p as 



N r fc-1 n 

*k— 1 



h=1 k=l,...,2N 



(21) 



As usual we rewrite the determinants in terms monic polynomials, and it follows directly 
from de Bruijn's integration theorem [3T] that the partition function (j!4p may be written as 
a Pfaffian, 

Z% m = (2N)\ Pf \\ [ d 2 z<(z)(,*-z)[^(z)^(z*)-^(,*)^(z)]l, (22) 

where p, ' (2) are monic polynomials of order k. It is clear that the expression (|22p for 
the partition function simplifies considerably if we choose the polynomials, p^,' m (z), to be 
skew-orthogonal with respect to the skew-product 

VWs = \f d 2 zw™(z)(z* - z) [f(z)g(z*) - f(z*)g(z)] . (23) 

Explicitly, the skew-orthogonality means that 

\P2k+l\P2i+l) S — \P2k \P2£ /S — U ' \ Zqa ) 

<P^i\P2t m )s = h k' m5 M, ( 24b ) 

where h 7 ^' 171 are constants. Using the skew-orthogonality of the polynomials we may write the 
partition function (|22p as |27| 

iV-1 

Z# m = (2iV)! Yl hl m . (25) 

fc=0 
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In the following we will always assume that the polynomials, p k ' m (z), are skew-orthogonal, i.e. 
that the polynomials satisfy the skew-orthogonality conditions (|24p . We will find the explicit 
realization of the polynomials, which are skew-orthogonal with respect to (|23p . at the end of 
this section. 

The skew-orthogonal polynomials are extremely useful, and not only the partition func- 
tion (|22p but also the eigenvalue correlations (|20p can be expressed in terms of the skew- 
orthogonal polynomials. One derivation of the eigenvalue correlations is given in |30| . which 
uses an idea presented in |32] , The eigenvalue correlations may be written as 



R% m (zi, ,z k ) = J] w™(z h )(zi - z h ) Pf 

h=l ~ J ~ 



n,m/ *\ n,rn/ * *\ 

K N \ z ii z j) K N \ z ji z i) 

n,mr \ n,m/ „, \ 

K N \ z ii z j) K N \ Z ji Z i) 



(26) 



where the so-called prekernel K^ m (u,v) may be expressed in terms of the skew-orthogonal 
polynomials, 

N—l n,m i \ n,m / \ n,m / \ n,m / \ 
.n.m,.. ^ ^P2k+l( U >P2k ( V )-P2k+l( V )P2k ( U 



K 



N 



I \ ■ p 2fe+lV u '^2fc \ u > i J 2k+l\ u )i J 2k V") / ~n 
(U,V)=2_^ JnTm • (27) 

fc=0 k 



In the literature one sometimes uses the notation of quaternion determinants, see e.g. [30] 123] . 
In this notation the eigenvalue correlations are written as 

k 

R n k ' m ( Zl , , z k ) = J] w™(z h )(z* h - z h ) Qdet K% m ( Zi , Zj ), (28) 
h=i 1<<J<* 

where Qdet is the quaternion determinant and K^ m (u,v) is a 2 x 2 matrix kernel. The 
relation between (|26p and (|28p should be obvious. This also explains why K^ m {u,v) is called 
the prekernel. The main point is that if we know the weight function, w™(z), and the 
prekernel, K^- m (u,v), then we know all the correlation functions. 

It only remains to determine the skew-orthogonal polynomials and the constants h^' m . 
For a general weight function the skew-orthogonal polynomials may be written as 

J&>) = (E[(*-^-<)) z n. m . ( 29a ) 

1=1 k 

=((*-x>- *)) ro* - *x* - z *)) x «, m > 

j=l i=l k 

which holds in complete analogue to [30j. It follows that the skew-orthogonal polynomials 
have real coefficients, Pk{z)* = pk(z*). In order to determine an explicit expression for the 
skew-orthogonal polynomials, we first notice that the weight function (|17p is invariant under 
rotation in the complex plane, hence w™{z) = if™(|z|). It is well-known that monomials are 
orthogonal with respect to such weights, i.e. we have the relation 

;>oo 

d 2 zw™(z)z k z* e = s n k ' m 5 M , with 4' m = 2vr/ dr <(r) r 2k+1 . (30) 

Jo 
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Using this relation it is straight forward to check that condition (|24ap is satisfied if 



P2kll(^ = aIld P2k m ( Z ) = ^Z C 2j 



m 2i 
Z J . 



(31) 



Furthermore, condition (124bp tells us that for any rotational invariant weight function, the 
constants c^j™ are given by a simple recursion relation, 



4j-l 



-2.i 



n,m f£°dr w™(r)r 
" 2j ~ 2 So 'drw™(r)r*j+ 1 ' 



(32) 



In the specific case where the weight function is given by equation (JTTJ) , the integrals can be 
performed analytically (see |29j ) and the skew-orthogonal polynomials are in monic normal- 
ization given by 



P2k+l( Z ) = Z 



2k+l 



k p k 

and p£f(*)=E II ( m + 2 iT 
1=0 \- j=t+\ 



■21: 



It is straight forward to check that with this choice of polynomials we have 



/ n,m I n,m 
\P2k+l\P2£ /S 



)S = h k °k£ 



■j-i i.n,m 

with n k 



i(7rr[m + 2A; + 2]) r 



(33) 



(34) 



Inserting the skew-orthogonal polynomials (|33p and the constants h^' m into equation (|27p we 
may write an explicit expression for the prekernel, 



(vrr[m + 1]) 



N-l k 
k=0 1=0 



,2k+l 



nj= (^ + 2j + 1)" K j=l {m + 2jY 



(u o v) 



■ (35) 



Note that the prekernel agrees with known results for n = 1 |3U] and n = 2 |23] , 

Due to formula (|26p we know all eigenvalue correlations, since we have explicit expressions 
for both the weight function (|17p and the prekernel 035 p . As an example, the density of states, 
or one-point correlation function, is given in terms of (|17p and (|35p as 



(36) 



Note that the first term on the right hand side of (|36p is responsible for a repulsion from the 
real axis. This repulsion can also be seen on figured] which shows scatter plots of the product 
matrix ([!]). 



4 Radial Behavior for the Density of States 

The radial behavior for the density of states for the product matrix, P™, is an important 
quantity for several reasons. This is true even though the density of states is not rotational 
invariant (see figure[T|). The radial behavior is important for quantities such as hole probability 
and overcrowding at the origin, since these quantities are related to the radial ordering of the 



8 




Figure 1: Scatter plots of 50 matrices defined as the product matrix, P™, with N = 25 
and m = 0, for n = 1 (left) and n = 3 (right). The eigenvalues have been rescaled by a 
factor (2N) n l 2 such that they approximately lie within the unit circle. Note that, due to 
the fact that eigenvalues come in complex conjugate pairs, the repulsion between individual 
eigenvalues results in a repulsion from the real axis. 



eigenvalues. Moreover, the radial behavior is often useful in comparison with numerical data. 
Therefore, in this section and section [5] we will mainly focus on the large- N behavior of the 
radial density of states. 

Note that studying the microscopic limits of the eigenvalue correlation functions, R 1 ^ ,m {z) 1 
in the complex plane for quaternion (/? = 4) Ginibre matrices is a challenging task; discussions 
of the microscopic origin limit have been given for n = 1 and n = 2 in [30j and [23j , respectively, 
but n > 3 remains an open problem. Averaging over the complex phase simplifies the formulae 
from the previous section considerably, which allows us to compare the results obtained in this 
paper directly to the results obtained for the product of complex (/3 = 2) Ginibre matrices 
in [E]. 

We define the radial density of states, Ppf m (r), as the average over the complex phase, i.e. 

P T(r) = ^fjORi' m (re i0 ), (37) 

where density of states, R™' m (z), is given by equation (|36p with weight function (jl 7[) and the 
prekernel fj35 1) . Due to the orthogonality relation (|30p it is straight forward to perform the 
integration over 6 in (|37p . We obtain 

2 / _ \ N ~ l r 4fc+2 

ffl-^U,..,™ r2 jE r[m + 2fc + 2] n - (38) 

k=0 

Note that the density of states is normalized to the number of eigenvalues rather than to 
unity, 

poo f 

dr2irrp^ m (r) = / d 2 zR^ m {z) = 2N. (39) 



9 



In the following we will discuss the large- N behavior of radial density of states. In the limit 
of large matrices, the eigenvalues are distributed within a circle with radius r ~ (m + 2N) n / 2 . 
We will first discuss the macroscopic density of states, then we will turn to the discussion of 
the microscopic behavior at the origin and at the edges as well as in the bulk of the spectrum. 

The explicit expression for the radial density of states (j38| for large TV" and r can be 
evaluated using the saddle point approach. A discussion of this is given in appendix [B] A 
similar discussion for the product of complex (/? = 2) Ginibre matrices was given in |15| . In 
the following we will also take m to be proportional to the size of the product matrix, i.e. we 



take m = 2Nm. For N » 1 and m n l 2 



< 



rn 



+ 2N) n l 2 we have (see appendix [B]) 



n,m i \ 

Pn ( r ) 




n r 2 / n - (m + 2N) 
2 rV™ 



erfc 



n r 2 l n 
2 



m 



l/n 



This structure leads us to introduce a scaled density of states defined as 



j% m (r) = (2A0"-VrV(2AT /2 ) 



(40) 



(41) 



We rescale of the radial variable, r, in order to get compact support, while the prefactor 
ensures that the density is finite in the limit of large matrices. Using the asymptotic behav- 
ior (|40p . the scaled density of states (|4ip can be written as 



Pn ( r ) 



T n 



±-2 



1 



mr 



erfc 



;2/n _ (tfr + 1) 



erfc 



m 



%1/n 



(42) 



The complementary error functions change only in the neighborhood of f = m n l 2 and f = 
(m + l)"/ 2 , respectively. Taylor expanding the arguments of the error functions around these 
values we obtain 



Pn ( r ) 




ANr- (m + 1)" /2 
~n (m + l)^- 1 )/ 2 



erfc 



4A^r 



m 



n/2 



(43) 



As a final remark in this section, we will determine the macroscopic (or mean) density of 
states. The macroscopic radial density of states is given by 



(r) 



km 

N^oo 



(r) 



-2. 
f n ' 



(e[(m + l) 



\n/2 



6[m 



n/2 



(44) 



where 9 denotes the Heaviside step function. Note that that macroscopic radial density of 
states for the product of n quaternion (J3 = 4) Ginibre matrices (|44p is similar to that complex 
(/3 = 2) Ginibre matrices, see |11[ [TSJ HU E]. The induced Ginibre ensemble has also been 
discussed for a single real (/3 = 1) matrix and a single complex (/3 = 2) matrix [23], which 
also gives a structure similar to (|44p with n = 1. The macroscopic radial density of states 044 p 
tells us that the eigenvalues of the product of induced Ginibre matrices lie within an annulus 
centered at origin in the complex plane. That this is an important generalization is illustrated 
by the single ring theorem |25| [26] . 



10 



0.3 



0.2 



0.1 









0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 

radial distance, r 



Figure 2: Plot for the radial density of states, p^} m (r), with N = 100 and m = 125 for n = 3. 
The dotted curve indicates the macroscopic limit. 

5 Microscopic Limits and Universality 

In this section we will discuss microscopic limits of the radial density of states. The asymptotic 
formulae of previous section enable us to identify the different regions where we can take the 
microscopic limit. We will discuss the universal bulk limit as well as the universality at the 
inner and outer edge, and finally we discuss the microscopic origin limit. 



Universal Bulk Limit The microscopic bulk limit can be taken, when the radial variable 
is far away from the edges, ttW 2 <C r <C (m + 2N) n l 2 . In order to study the microscopic 
behavior we must rescale the variable r — > r(2N)~ n / 2 to get compact support, and then make 
another rescaling r — > r(2N) n l 2 to get to the level of the eigenvalue fluctuations. Clearly, the 
two scalings cancel each other. The final step to obtain the microscopic bulk limit is to unfold 
the variable, such that the density of states becomes locally flat. From previous section we 



know that the unfolding can be made with a change of variable given by f 
region of interest, m n l 2 <C r <C (m + 2N) n / 2 , this means that 



- 2 / n . In the 



dr = 2irr 



v n,m/ \ 

hm prf (r) 



N-s-oo 



dr, 



which indeed gives a flat density. We get a microscopic bulk radial density given by 



n,m / ~\ 



i 



(45) 



(46) 



Note that the microscopic bulk limit of the radial density of states is independent of n and 
m. Moreover, the limit is also equivalent to the microscopic bulk limit of the radial density of 
states for the product of complex (/? = 2) Ginibre matrices |15| . Of course the main interest 
in the bulk is the higher order correlation functions in the complex plane, R?' m (z). 
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Universal Edge Limit We will now investigate the microscopic behavior at the inner and 
outer edge, i.e. for r ~ m™/ 2 and r ~ (m + 2N) n / 2 , respectively. Looking at equation (j43]) 
from previous section, we introduce the edge variables, £\ n and £? out , given by 

2Nf-m n / 2 , [2N r - Cm + lW 2 , . 

and e out = \ — , ■ , u „_ u/9 , (47) 



n m (n-i)/2 OUI V n (m + l)^- 1 )/ 2 

with the scaled quantities f and m defined in previous section. The error functions changes 
only in a narrow region, and the width of this region is proportional to l/y/2N when we use 
the scaled variable f. With this in mind, it clear that the edge variables (|47p zoom in on the 
region of change (for the inner or outer edge, respectively) and then rescale such that the 
edge is located at unity. 

With the above definition of the edge variables (j47| we get the universal edge behavior 

) = J im Pn" 1 ^ = T- ^[V2i], (48) 

where e is defined as i\ n and i out for the inner and outer edge, respectively. The result holds 
for any n and rh and is thus universal; see also the discussions given in |27| |22| 124"! 115) . It 
should be mentioned that if m is small, i.e. m is not proportional to the size of the matrix, 
then there is no inner edge. Instead there exists a microscopic origin limit which does depend 
on both n and m; we will discuss this limit below. 

Microscopic Origin Limit We have already seen that when m ^> 1, then the eigenvalues 
are repulsed from the origin. Here we will look at the situation where m is of order unity, 
and discuss the microscopic behavior at the origin. The microscopic limit can be obtained by 
taking N — > oo while keeping the radial variable far away from the edge, r <C {2N) n l 2 . If we 
do so directly in equation fj38|) we get 

£SnM - J™ PT(T) = - G n ; n [ m ^^ m r 2 ) g r[m + 2fc + 2]" - (49) 

The infinite series can be recognized as a generalized hypergeometric function, and we can 
write the microscopic origin limit as 



n,m / \ ^' /~in,0 

PoriginlH - nT [ m + 2 ] 0,ny m ^ ^ m 



x 1^2n( m+2 m+2 m+3 m+3 



4 



^ J . (50) 

where both (m + 2) /2 and (m + 3) /2 appear n times in the hypergeometric function. Note 
that the microscopic origin limit depends on both n and m. This is to be compared with the 
microscopic limit for the product of complex (/3 = 4) Ginibre matrices given in |15] , 
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6 Conclusions and Outlook 



In this paper we have investigated the spectral properties for the product of n independent 
induced quaternion (/? = 4) Ginibre matrices. We have shown that the joint probability 
density function for the product matrix is identical to that for a single Ginibre matrix, except 
that the weight function is more complicated. Furthermore, we have seen that the weight 
function for the product quaternion (/3 = 4) Ginibre matrices is identical to that for the 
product of complex (/3 = 2) Ginibre matrices, i.e. is given in terms of a Meijer G-function. 

For matrices of finite size, we have used the method of skew-orthogonal polynomials to 
obtain an explicit expression for the prekernel. This implies that we can write down any 
correlation function of a product matrix with finite size, since all eigenvalue correlations can 
be expressed in terms of weight function and the prekernel. 

Finally, we looked at the radial density of states, defined from the density of states in 
the complex plane by averaging over the complex phase. In this case we investigated the 
limit of large matrices. We found that the macroscopic radial density was similar to the 
radial behavior for the product of complex (/? = 2) Ginibre matrices. Furthermore, we 
discussed the microscopic limits. For the weakly induced ensemble, m = O(N ), we studied 
the microscopic origin limit, while for the strongly induced ensemble, m cx N, we discussed 
the universal behavior at the inner and outer edge. We also looked at the behavior in the 
bulk. We saw that the edge behavior, given in term of a complementary error function, was 
universal and identical to the case of complex Ginibre matrices. The microscopic origin limit 
was written using a generalized hypergeometric function. 

It would be a natural extension of this paper to study the large matrix limit of the 
correlation functions in the complex plane, i.e. without taking the average over the complex 
phase. Furthermore, it is an appealing challenge to extend the discussion of products of 
Ginibre matrices to include the case of real (/3 = 1) Ginibre matrices. We expect that the 
results for the edge and the bulk to hold for the complex eigenvalues, while the origin will be 
generalized as already known for the product of two matrices. It would also be interesting to 
discuss the hole probability and overcrowding at the origin for quaternion matrices, this was 
done for the product of complex matrices in |16j . Diagrammatic methods have previously been 
used to discuss the macroscopic limit for the product of complex rectangular matrices |12| 
and it would be desirable to understand this on the level of orthogonal polynomials. Several 
of these projects are currently on the way. 

Acknowledgment G. Akemann, M. Kieburg, E. Strahov as well as the participants at 
the 'VIII Brunei-Bielefeld Workshop on Random Matrix Theory' are thanked for useful dis- 
cussions. The author is supported by the German Science Foundation (DFG) through the 
International Graduate College 'Stochastics and Real World Models' (IRTG 1132) at Bielefeld 
University. 

A Derivation of the Jacobian 

In this appendix we calculate the Jacobian for the change of variables corresponding to the 
parameterization given in section (|2|), i.e. the the change of variables from to Aj, Tj and 
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Recall that X, are quaternion (/? = 4) Ginibre matrices, Aj are diagonal (block) matrices, Tj 
are strictly upper-triangular (block) matrices and £7j 6 Sp(N)/U(l) N are symplectic matrices. 

In the following dX = [dX ab ] denotes one-forms ordered as a matrix, while DX = f\ dX ab 
denotes the corresponding exterior product of independent one-forms, and \DX\ is the un- 
orientated volume element. The one-forms of dXi are related to those of dAj, cffj and dlli 
as 

dXi = UidYJjr^, (51) 

where 

with dMi defined as 

{dMi) ab = (dA t ) ab {Ai) bb - {Ki) aa {dA i+1 ) ab + ^{dAi) ac {Ti) cb - ^(T i ) ac (cL4 i+1 ) c6 . (53) 



dYi = dAi + dTi + dM u (52) 



N N 



c=l c=l 



Here dA n+ \ = dA\ and tL4j = U^ 1 dUi are anti-self-dual matrices [27]. Due to the anti- 
self-duality the one-forms in the upper and lower triangular entries are related by (dAi) ab = 
— (dAi) ba , where the bar denotes the quaternion conjugate |27) . For this reason we will only 
keep the lower triangle as independent one-forms. Furthermore, it should be noted that since 
U{ S Sp(N)/U(l) N the diagonal entries (dAi) aa has two independent one-forms (each entry 
in the lower triangle has four independent one-forms), 



U/.l,)„„: mPU + HdAW) aa = ±( d ° ia d <). (->■.]) 



where a ia = -(A^ ) aa + i(A^) aa G C. 

We can now calculate the Jacobian. It is trivially seen that the Jacobian for the change 
of variables from X to Y equals one, since Ui are symplectic matrices. In order to find the 
Jacobian for the change variables from Y to A, T and M we write (|52[) as 

(dY} q) ) ab = (dA® ) ab + {dM^U for a = 6 and g = 0,l, 

{dY % {q) ) ab = (dT^U +{dMl q) ) ab for a< 6 and g = 0,l, 

(dY l (q) ) ab = (dT l {q) ) ab + {dM^) ab for a< 6 and g = 2,3, 

(dY} q) ) ab = + ( dM ( q) ) ab for a > b and 9 = 0, 1, 

(dY^ q) ) ab = (dM^U for a> b and 9 = 2,3. (55) 

With this ordering it is straight forward to see that the Jacobian matrix for the change of 
variables Y to A, T and M is upper triangular and that the Jacobian is equal to unity, such 
that 

\DX\ = \DY\ = \DA\\DT\\DM\. (56) 

Here DA and DT are the Euclidean volume elements for A and T, respectively. The measure 
DM is given by 

TV 3 n N 1 n 

DM = ( A A A ( dM i q) U) A (A A MdM^) bb ). (57) 

a>bq=0i=l fe=lq=0i=l 
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Note that the volume elements DM and DA contain the same number of independent one- 
forms, which was to be expected from (|53p . As the final step we want change variables from 
M to A. First we will show that the Jacobian matrix J = dM/dA can be chosen to be 
triangular and that the Jacobian \J\ is independent of T. To do so, we order the matrix index 
a = ab such that ab < a'b' if a > a' or if a = a' and b < b'. This is identical to the ordering 
used in [15] . The one-forms dM and dA are related by 

(dM i ) a = J2( J ij)aa'(dA j ) a ,, (58) 

where J is the Jacobian matrix. Using that the matrices Tj are upper triangular, we see 
from (|53p that the entries with a > a' in the Jacobian matrix are zero, i.e. the Jacobian 
matrix is triangular such that only the diagonal entries (a = a') contribute to the Jacobian, 
\J\. It is clear from fj53|) that the diagonal entries do not depend on T, and we can therefore 
replace the DM with 

(dMi) ab = (dAi) ab (Ai) bb - {A t ) aa {dA i+l ) ab . (59) 

Using (|59p in (|57p it is straight forward to obtain the Jacobian. Recall that the entries (dMi) ab 
and {dAi) ab are quaternions, or equivalently 2x2 matrices. For a > b we have that 

\d Vl du*J ab \0 xj bb \0 xj aa \dv i+1 du i+1 ) ab 

(60) 










* 


Xiadu^ 


-l)ab 




l)ab 



XibdUiab Xia.du^+i'jai) x i b d v i ab + X io.dv*^ + -^ ab 

Xibdv iab - x* ia dv (i+1)ab x* ib du* ab - x* a du* {i+1)ab 



A similar equation holds for a = b but in that case there are no diagonal terms. Instead of 
using the quaternion index (q), we will use the matrix index k£, i.e. (dMj)^ with k,£ = 1,2 
and similarly for dA. We see that 

n n 

(a>b) : f\ {dMi)H = (x lb x 2b ■■■x nb - x la x 2a ■ ■ ■ x na ) f\ du iab , 
i=i i=i 

n n 

(a>b) : f\ (dMi)H = (x lb x 2b ■■■x nb - x\ a x\ a ■ ■ ■ x* na ) f\ dv iab , etc. (61) 
i=i i=i 

Combining these results we find the Jacobian 

N N 



\J(z)\ 



dM(z) 



dA 



Y[ \ Z b ~ Za\ 2 \zb ~ Z*\ 2 Y[ W ~ Z*\ 2 , (62) 
a>b c=l 



where z is a short hand notation for the product z a = Y\7=\ x i a - This is the Jacobian given 
in section [2] 
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B Asymptotic of Radial Density of States 



In this appendix we discuss the asymptotic behavior of the radial density of states. We study 
the asymptotics for large matrices (N 3> 1) and large r (this means 1 <C r J$ (2N) n l 2 if m 
is of order unity and m™/ 2 J$ r J$ (m + 2N) n l 2 if m > 1). Recall that the radial density of 
states is given by 



" it ' V m, . . . , m 



,2 



N-l 



4k+2 



^ r[m + 2fc + 2]"' ^ 



The asymptotics for the Meijer G-function at large arguments is known, see [33]. Furthermore, 
the particular Meijer G-function given in (|63p has been discussed explicitly in [15] . For this 
reason we refer the reader to the above mentioned references and state the result without 
derivation, 



,^n, o 

^0 n ■ 



r 2 ) « (27F)( ^ l m r d-n)/n r 2m e -nr^ for r ^ \, ( 6 4) 

n 



In order to evaluate the sum in (|63p we will use a saddle point approximation. For large N the 
sum can be replaced by an integral. If we also make a shift in the integration, k — > (k—m—1) /2, 
we see that 

N ~ 1 ^4fc+2 i r m+2N J2k-2rn 

V ~ - / rfjfc J! (65) 

^ T\m + 2k + 2} n 2 L IU + ll n V 7 



fc=0 



The gamma function is well approximated by Stirling's formula, T[A;+1] ~ v / 2irk{k / 'e) k '. After 
this replacement we have 

iY " 1 r 4k+2 (n \-n/2 f m+2N 

y ^ r , r « ^ r ~ 2m / ^ A;""/ 2 e 5 ^ , (66) 



k=0 

where 



S(k) = k log r 2 — n£; log k + n/c. (67) 



Note that a similar integral also appeared in the evaluation of the kernel for the product of 
complex {(3 = 2) Ginibre matrices [15] . In order to evaluate the integral we approximate the 
exponential in (j66]) by a Gaussian function. The solution of the saddle point equation, 

S'(K) = log r 2 -n log h = 0, (68) 

gives the location of the maximum of the Gaussian function, k* = r 2//n . At the maximum we 
have S{k*) = nk* and S"(k*) = —n/k*, and therefore we have that 

Y , r , - « W r ~2m e nk, / dk k -n/2 n{ k-K? /2K (6g) 

^r[m + 2A; + 2]« 2 J m y ' 
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We have r> 1 such that the integrand is a narrow peak centered at k*. For this reason we 
can replace the integration limits with ±00 if m <C fc* <C m + 2N . In this region it straight 
forward to perform the integral, 



jV-l 

E 

fc=0 



^ik+2 



( 2vr )(i-«)/2 



r[m + 2A; + 2] n 



r -2m r (l-n)/n e nr 2 /" 



(70) 



Inserting this result and the asymptotic behavior of the Meijer G-function ()64[) into (|63p . we 
obtain 



Pn ( r ) 



1-2 



for m n/2 < r <C (m + 2iV) 



n/2 



(71) 



This gives the large- iV asymptotic of the radial density of states in the bulk, but if we want to 
study the behavior at the edges we need to be more careful. If the maximum of the Gaussian 
peak is in the neighborhood of the edges of the integration interval, ~ m and ~ m + 2N, 
then we cannot replace the integration limits with ±00. Instead we must consider situation 
where the Gaussian peak lies partly outside the integration interval. If we write (j69]) as 



N-l 

E 

k=0 



„4fc+2 



(2*) 



-n/2 



T[m + 2k + 2] n 



r -2m e nk* 



m+2N 



dk 



dk\k: n/2 e- n{ - k - k ^' 2k % 



(72) 

it is clear that the evaluation of the integral at the inner and outer edge is identical, since 
the Gaussian peak is symmetric around its maximum. Partly integrating over the Gaussian 
peak results in a complementary error function. Combining this with the discussion above 
we obtain the following asymptotic behavior of the radial density of states, 



n,m 1 \ 

Pn ( r ) 



f n 



1 



nir 2 



- erfc 



r 2/« _ ( m + 2N) 



,1/n 



erfc 



n r 2 l n — in 
2 



,1/n 



(73) 



If m is small (m ~ 1 <C N), then we only have an outer edge and the density of states becomes 



Pn ( r ) 



T n 



1 



nir 2 



erfc 



n r 2 ' n - 2N 



,1/n 



(74) 



which holds for 1 <r ^ 2N . 
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